Signal to Noise Ratio Determination
Difficulty Level:
Tags pre-process☁quality☁noise☁snr

SNR, standing for Signal to Noise Ratio , is an estimation of how much noise there is in a digital signal. All digital signals can be divided in two main elements: the relevant data that we meant to record and the non relevant noise that got unavoidably recorded along with it. An ideal signal would have a very high SNR value. The higher the value, the less noise there is in the signal. However this is not always easy to achieve in real life.

Because the presence of noise usually shadows the "real" and relevant data we want to analyse, it is important to be able to filter as much noise as possible without unwillingly removing relevant pieces of information. There is already a really interesting Jupyter Notebook explaining the process of digital filtering electrophysiological signals, which you can find here .

In this Jupyter Notebook , however, we will focus on how to estimate the SNR in a digital signal, therefore determining the quality of the recording.

1 - Importation of the needed packages

In [1]:
# biosignalsnotebooks own package for loading and plotting the acquired data
import biosignalsnotebooks as bsnb

# Scientific packages
from numpy import array, linspace, mean, log10
from scipy.signal import butter, filtfilt
In [2]:
# Base packages used in OpenSignals Tools Notebooks for plotting data
from bokeh.plotting import figure, output_file, show, curdoc
from bokeh.io import output_notebook
from bokeh.layouts import column
from bokeh.models import ColumnDataSource, Plot, LinearAxis, BoxAnnotation, Arrow, VeeHead, LinearAxis, Range1d
output_notebook(hide_banner=True)

2 - Load of sample signal data

In [3]:
# Load of data (using a relative path to the project folder)
data, header = bsnb.load("../../signal_samples/eeg_sample_artefacts_seg1.h5", get_header=True)

3 - Storage of sampling frequency and acquired data inside variables

Since in this case we know that the original signal is an EEG recording, we can also convert its units to $\mu V$ following the method explained in the EEG Sensor - Unit Conversion notebook.

In [4]:
# Sampling frequency and acquired data
fs = header["sampling rate"]

# Let's get de raw signal data and the time range of the recording
channel = list(data.keys())[0]
signal_raw = data[channel]
time = linspace(0, len(signal_raw) / fs, len(signal_raw))

# Let's convert the signal's units, since we know it is a ECG signal
vcc = 3e6 # uV
gain = 40000
resolution = header['resolution'] # Resolution (number of available bits)
signal = (((array(signal_raw) / 2**resolution) - 0.5) * vcc) / gain
In [5]:
# Let's plot our raw signal
p = figure(plot_width=1200, plot_height=200, x_axis_label="Time (s)", y_axis_label="EEG (\u03BCV)", **bsnb.opensignals_kwargs("figure"))
p.ygrid.grid_line_alpha=0.5

# add a circle renderer with x and y coordinates, size, color, and alpha
p.line(time, signal, **bsnb.opensignals_kwargs("line"))

# apply OpenSignals style
bsnb.opensignals_style([p])

show(p) # show the results

4 - SNR estimation

This notebooks" proposed approach to estimate the SNR consists of four phases:

  • Smoothing the original signal to remove as much noise as possible but retaining our signal component;
  • Subtracting the smoothed signal from the original unfiltered acquisition, which isolates our noise component;
  • Determining the signal and noise peak-to-peak amplitude;
  • Estimating the SNR value by dividing the previously determined peak-to-peak amplitude.

4.1 - Signal smoothing

This step depends much on the characteristics of the signal we are dealing with. Different types of signals and acquisitions with varying degrees of quality often require different strategies to effectively reduce the noise components without losing relevant information. In the same way, the quality and accuracy of the filtering will determine the SNR value.

In this tutorial we will follow a conservative EEG filtering using a Butterworth bandpass filter to remove frequencies lower than 1 Hz and higher than 50 Hz, which are usually outside the range of brainwave frequencies. There is, however, a more detailed Jupyter Notebook dealing with the details of EEG filtering: Digital Filtering - EEG .

In [6]:
# Time window to center our processing
sample_start = 0
sample_end = 30*fs
# Baseline shift of window
signal_shift_window = array(signal[sample_start:sample_end]) - mean(array(signal[sample_start:sample_end]))

# Butterworth bandpass filter
low_cuttoff_wide = 1# lower cutoff frequency for bandpass filter (Hz)
high_cuttoff_wide = 50 # Upper cutoff frequency for bandpass filter (Hz)

# Creating a bandpass Butterworth filter
b, a = butter(2, [low_cuttoff_wide, high_cuttoff_wide], btype='bandpass', fs=fs)
bandpass_filtered = filtfilt(b, a, signal_shift_window) # filter with filtfilt to obtain a zero phase response 
In [7]:
# Let's plot our raw signal
p = figure(plot_width=1200, plot_height=200, x_axis_label="Time (s)", y_axis_label="EEG (\u03BCV)", **bsnb.opensignals_kwargs("figure"))
p.ygrid.grid_line_alpha=0.5
p.title.text = 'Comparison of the original and filtered signals'

# add a circle renderer with x and y coordinates, size, color, and alpha
p.line(time[sample_start:sample_end], signal_shift_window, line_alpha=.7, legend_label='Original signal', **bsnb.opensignals_kwargs("line"))
p.line(time[sample_start:sample_end], bandpass_filtered, legend_label='Filtfilt bandpass filter', **bsnb.opensignals_kwargs("line"))

p.x_range = Range1d(-5, time[sample_end+7000])

# apply OpenSignals style
bsnb.opensignals_style([p])

show(p) # show the results

# Closer look to the efects of the filter
p.title.text = 'Closer look'
p.x_range = Range1d(time[9000], time[11500])
p.y_range = Range1d(-10, 15)

# apply OpenSignals style
bsnb.opensignals_style([p])

show(p)

4.2 - Subtracting filtered signal to the original signal

In order to isolate the noise removed in the filtering process.

In [8]:
# Let's subtract our filtered signal from the original signal 
noise_component = signal_shift_window - bandpass_filtered
In [9]:
# Let's plot our raw signal
p = figure(plot_width=1200, plot_height=250, x_axis_label="Time (s)", y_axis_label="EEG (\u03BCV)", **bsnb.opensignals_kwargs("figure"))
p.ygrid.grid_line_alpha=0.5
p.title.text = 'Comparison between the original signal, the filtered signal and the subtracted noise'

# add a circle renderer with x and y coordinates, size, color, and alpha
p.line(time[sample_start:sample_end], signal_shift_window, line_alpha=.5, legend_label='Original signal', **bsnb.opensignals_kwargs("line"))
#p.line(time[sample_start:sample_end], filtered_signal_wide, line_color='red', legend_label='Wide passband filter')
p.line(time[sample_start:sample_end], bandpass_filtered, line_alpha=.5, legend_label='Filtfilt bandpass filter', **bsnb.opensignals_kwargs("line"))
p.line(time[sample_start:sample_end], noise_component, legend_label='Noise component', **bsnb.opensignals_kwargs("line"))

p.x_range = Range1d(-5, time[sample_end+7000])

# apply OpenSignals style
bsnb.opensignals_style([p])

show(p) # show the results

# Closer look to the effects of the filter
p.title.text = 'Closer look'
p.x_range = Range1d(time[2000], time[3500])
p.y_range = Range1d(-10, 15)

# apply OpenSignals style
bsnb.opensignals_style([p])

show(p)

The big oscillations we can see in the noise component correspond to the lowest filtered frequencies (< 1 Hz). If you take a closer look at the red signal, you can see that it actually contains many small high frequency oscillations too, which correspond to the high filtered frequencies (> 50 Hz).

4.3 - Calculating the signal"s and noise"s peak-to-peak amplitudes

In [10]:
ptp_signal = max(signal_shift_window) - min(signal_shift_window)
ptp_noise = max(noise_component) - min(noise_component)
In [11]:
print('Signal\'s component peak-to-peak amplitude:', ptp_signal, '\nNoise component\'s peak-to-peak amplitude:', ptp_noise)
Signal's component peak-to-peak amplitude: 22.5677490234375 
Noise component's peak-to-peak amplitude: 7.65416513529668

4.4 - Estimating the SNR value using the peak-to-peak amplitudes

The calculation of the SNR is as follows: \begin{align} \mathrm{SNR_{dB}} = 10 \log_{10} \left ( \frac{P_{signal}}{P_{noise}} \right ).\\ \end{align}

where $P_{signal}$ stands for the power of the original signal and $P_{noise}$ for the power of the noise component.

The logarithm allows to convert the SNR ratio, which has no units, into decibels (dB). This is done for convenience, because the SNR can have very large values. Converting the ratio to a logarithmic scale makes it much easier to manage.

In [12]:
# Estimation of SNR using the peak-to-peak amplitude of the signal and the noise component.
snr_peak_to_peak = 10*log10(ptp_signal/ptp_noise)
In [13]:
print('SNR(dB): ', snr_peak_to_peak)
SNR(dB):  4.695904156683456
In [14]:
from numpy import var, sqrt
# HIDDEN: Alternative SNR estimation methods
# Estimation of SNR using the variances of signal and noise component
snr_variance = 10*log10(var(signal_shift_window) / var(noise_component))

#print('SNR(dB): ', snr_variance)

# Estimation of SNR using the RMS of signal and noise component
rms_signal = sqrt(mean(signal_shift_window**2))
rms_noise = sqrt(mean(noise_component**2))
snr_rms = 10*log10((rms_signal/rms_noise)**2)

#print('SNR(dB): ', snr_rms)

A SNR value higher than 0 dB indicates that there is more relevant signal than noise. However, in this case, although the value is positive, it is also low. That means that the filtered noise is significantly shadowing the real information. Nonetheless, in this case the SNR value will change depending on the cutoff frequencies chosen for the bandpass filter, so it is very important to analyse our data and chose the optimal way to extract the noise. Obviously, more robust filtering methods ensure a more accurate estimation.

Furthermore, it is necessary to differentiate between filtering used to remove noise and filtering used to extract the data we want to focus in. In the case of EEG, for example, the analysis of the alpha bands usually requires the isolation of frequencies between 8 and 13 Hz, but that does not mean that everything outside that range is noise. There may be other real information that is not useful for our current analysis and therefore must be ignored, but not classified as noise. So it is always important to determine what is to be considered as noise in the signal we are working with.

In this example, the cutoff frequencies where chosen because in EEG, frequencies below 1 and higher than 50 Hz are prone to contain more noise than useful data (with exceptions, of course).

Usually the acceptable value for SNR depends on the type of the signal and the use we are going to give to it, but as a general rule of thumb the higher it is, the better the quality of signal.

We hope you have enjoyed this tutorial!
You can keep learning thanks to the other available Jupyter Notebooks that explain how to process, analyse and extract features from electrophysiological signals.
Check the list of all available notebooks here

In [15]:
bsnb.css_style_apply()
.................... CSS Style Applied to Jupyter Notebook .........................
Out[15]: